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1.0 

INTRODUCTION 

This  report  summarizes  work  done  on  contract  N00014-73-C-0458 
during  the  period  3-1-78  to  4-30-80,  concentrating  especially  on  the 
parts  of  the  program  dealing  with  the  marine  environment  during  the  last 
year  of  the  contract.  Efforts  on  the  beach  environment  tasks  during 
the  last  year  are  summarized  in  a  separate  report  (no.  134400-12-F) , 
and  results  obtained  during  the  first  year  of  the  contract  are  described 
in  report  no.  134400-7-T.  Various  aspects  of  the  technical  work  have 
also  been  reported  in  papers  submitted  to  Remote  Sensing  of  Environment 
and  the  IEEE  Transactions  on  Geoscience  and  Remote  Sensing.  The  tech¬ 
nical  monitor  for  this  contract  was  Hr.  Hans  Dolezalek.  The  principal 
investigator  was  Dr.  David  Lyzenga  and  the  project  manager  was  Mr.  Fred 
Thomson. 

Three  areas  of  research  are  described  in  this  report.  The  first 
deals  with  the  extraction  of  water  optical  properties  from  remote 
sensing  data,  including  both  passive  (satellite  and  aircraft)  data  and 
active  (aircraft)  data.  These  properties  are  needed  for  use  by  the 
bottom  recognition  and  water  depth  algorithms,  and  may  have  independent 
applications  as  well.  The  second  topic  is  a  further  evaluation  of  the 
bottom  recognition  algorithm  using  both  Landsat  and  aircraft  data  over 
a  test  site  near  North  Cat  Cay  in  the  Bahamas.  The  third  topic  deals 
with  radiative  transfer  modeling  as  applied  to  water  and  atmospheric 
problems.  A  Monte  Carlo  model  has  been  developed  which  is  capable  of 
accurate  solutions  and  of  modification  for  complex  geometries  and 
boundary  conditions,  but  requires  a  fairly  large  expenditure  of  compu¬ 
ter  time.  Secondly,  a  re-evaluation  of  several  approximate  models  has 
been  conducted,  including  the  quasi-single-scattering  (QSS)  model  which 
has  bBen  used  in  the  past  for  water  reflectance  calculations  and  the 
double-delta  (Turner)  model  which  has  been  used  for  atmospheric 
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calculations.  Recently  a  rather  serious  deficiency  has  been  discovered 
in  the  Turner  model  which  causes  an  overestimate  of  the  path  radiance 
in  the  antisolar  direction.  The  QSS  model  appears  to  give  adequate 
results  for  the  conservative  scattering  case  with  small  optical  thick¬ 
ness  but  breaks  down  for  larger  values  of  the  optical  thickness.  A  new 
model  has  been  developed  which  shares  certain  similarities  with  both 
the  Turner  and  QSS  models  but  gives  somewhat  better  results  for  the 
conservative  scattering  case,  particularly  for  large  optical  depths. 
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2.0 

WATER  PROPERTIES 

2.1.  PASSIVE  METHODS 

The  basis  for  the  methods  of  processing  passive  shallow-water 
remote  sensing  data  described  in  this  report  is  the  assumption  of  an 
exponential  depth  dependence  of  the  observed  bottom- reflected  radiances. 
The  rationale  for  this  assumption  has  been  discussed  in  previous 
reports  [1,2],  and  no  overwhelming  evidence  has  been  found  to  the  con¬ 
trary,  although  it  is  certain  to  break  down  in  areas  where  the  optical 
properties  of  the  water  are  changing  rapidly.  In  areas  where 
horizontal  mixing  ensures  a  relatively  homogeneous  composition  of  the 
water,  however,  the  assumption  appears  to  be  justified.  Thus,  when  the 
measured  radiances  are  transformed  according  to  the  following  equation 

X.  =  ln(L.  -  L  .)  (1) 

1  1  si 

where  L,  is  the  measured  radiance  in  band  i  and  L  .  is  the  deep-water 
1  si 

radiance,  the  resulting  variables  are  linear  functions  of  the  water 

depth  and  are  linearly  related  to  each  other,  for  a  given  bottom 

reflectance.  That  is,  if  is  plotted  versus  X^  and  the  water  depth 

is  varied,  the  data  points  will  fall  along  a  straight  line  whose  slope 

is  K  /K  ,  where  K  is  the  irradiance  attenuation  coefficient  of  the 
i  J  i 

water  in  band  i.  If  the  bottom  reflectance  is  changed,  the  data  points 
will  fall  along  a  parallel  line  which  is  displaced  from  the  first  by 
an  amount  proportional  to  the  change  in  reflectance. 

This  suggests  a  procedure  for  obtaining  water  attenuation  para¬ 
meters  from  shall cw-water  MSS  data.  If  a  set  of  radiance  measurements 
are  made  over  an  area  of  variable  depth  but  uniform  bottom  reflectance, 
a  linear  correlation  would  be  expected  between  any  two  of  the  variables 
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X  and  X  .  A  regression  analysis  of  this  data  would  then  yield  the 
>  .1 

ratio  of  the  attenuation  coefficients  in  the  two  wavelength  bands  used. 
This  procedure  does  not  yield  the  correct  results,  however,  if  the 
bottom  reflectance  varies  with  depth  within  the  area  selected.  Thus, 
the  problem  is  to  select  an  area  of  uniform  bottom  reflectance  for  the 
analysis.  This  selection  process  is  directly  analogous  to  the  process 
of  training  set  selection  in  conventional  supervised  classification 
techniques  used  in  remote  sensing.  The  selection  may  be  done  on  the 
basis  of  field  observations  and/or  photographic  interpretation  using 
color  aerial  photography.  Alternatively,  an  analysis  of  the  self- 
consistency  of  the  data  can  be  done  to  determine  the  uniformity  of  the 
training  set.  In  this  case,  a  uniform  bottom  is  indicated  by  a  high 
correlation  among  the  X.  variables,  and  this  condition  can  usually  be 
identified  by  an  inspection  of  scatter  plots  of  the  data.  Figure  la 
shows  a  scatter  plot  of  values  over  an  area  containing  a  mixture  of 
bottom  types,  while  Figure  lb  shows  the  distribution  of  the  same 
variables  over  a  uniform  sand  bottom. 

Having  selected  a  training  set,  a  conventional  least-squares 
regression  analysis  may  be  done  to  determine  the  ratio  of  attenuation 
coefficients.  The  disadvantage  of  this  ;  rocedure,  however,  is  that 
the  results  depend  upon  the  order  in  which  variables  are  selected, 
since  the  mean  square  deviation  from  the  regression  line  is  calculated 
in  the  direct  ion  of  the  dependent  variable.  A  better  procedure  is  to 
determine  the  line  for  which  the  mean  square  deviation  measured  perpen¬ 
dicular  to  the  line  is  minimized.  This  procedure  results  in  the 
equat ion 

K  /K  =  a  +  V  +  1  (2) 

i  .1 

who  re 

a  -  41— "JJ-  (3) 

>  J 


A 


and  o..  =  X.  X  -  X  X  (4) 

1 J  i  .1  i  .! 

(.  i . o .  ,  '  .  .  is  tlie  variance  of  the  X.  measurements,  o  is  the  variance 

1  jj 

of  the  X.  measurements,  and  n . .  is  the  covariance  of  X  and  X  ).  This 
J  t.l  t  j 

equation  yields  a  result  which  does  not  depend  on  the  order  in  which 

the  variables  are  selected. 

This  procedure  has  been  tested  with  two  data  sets  collected  near 
North  Cat  Cay  in  the  Bahamas.  The  test  site  selected  for  this  evalua¬ 
tion  is  located  on  the  western  edge  of  the  Great  Bahama  Bank,  abqnt 
90  km  east  of  Miami,  Florida,  and  15  km  south  of  Bimini.  The  bottom  in 
this  area  is  composed  largely  of  white  carbonate  sand  with  some  rock 
and  coral  formations,  and  heavy  growths  of  vegetation  (mainly  Thai  ass  la) 
in  protected  areas.  The  water  is  very  clear,  and  extensive  shallow 
areas  exist  with  depths  ranging  from  1  to  15  meters,  dropping  off  to 
very  great  depth  at  the  edge  of  the  Batik.  A  map  of  the  test  site  is 
shown  in  Figure  2. 

Measurements  of  the  irradiance  attenuation  coefficient  of  the 
water  were  made  at  several  locations  in  the  Great  Bahama  Bank  in 
October  1977,  using  Kahl  submarine  photometer  with  filters  approxi¬ 
mating  the  Landsat  green  (MSS4,  0.5-0. 6  pm)  and  red  (MSS5,  0.6-0. 7  pm) 
bands.  Two  of  these  measurements  were  taken  near  the  North  Cat  Cay 
test  site,  on  October  4  and  5,  1977.  Landsat  data  was  also  collected 
over  the  test  site  on  October  11,  1977.  The  Landsat  MSS  data  was 
collected  in  high  gain  mode  and  exhibits  a  large  amount  of  striping, 
which  was  corrected  prior  to  the  processing  described  in  this  report. 

The  test  site  was  revisited  in  August,  1978,  at  which  time  a  set 
of  observations  of  the  bottom  type  and  water  depth  were  made  at  nine 
stations  around  North  Cat  Cay.  Underwater  photographs  were  taken  and 
were  subsequently  analyzed  in  order  to  extract  the  bottom  reflectance 
at  each  station.  Aircraft  MSS  data  was  also  collected  at  tills  time 
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using  the  ERIM  H-8  mul t ispect ra 1  scanner  system.  This  data  was 
collected  at  an  altitude  of  3.3  km,  with  a  spatial  resolution  of 
approximately  8  meters,  as  compared  to  the  Landsat  spatial  resolution 
of  79  meters.  The  spectral  resolution  of  the  aircraft  scanner  was  also 
somewhat  higher  than  the  Landsat  MSS,  with  bandwidths  on  the  order  of 
0.05  nm  for  the  wavelengths  bands  used.  The  aircraft  data  included 
some  sun  glint  and  path  radiance  effects,  which  were  removed  during  a 
pre-processing  step. 

The  aircraft  MSS  records  a  digital  signal  which  is  proportional  to 
the  radiance  at  each  point  in  the  scene.  The  radiance  within  each  of 
the  spectral  bands  indicated  in  Table  1  is  recorded  separately,  and 
with  a  different  proportionality  constant  depending  on  the  gain 
settings  selected  by  the  operator.  The  proportionality  constants 
relating  the  digital  signals  to  the  radiances  in  each  band  could  be 
determined  by  an  analysis  of  the  calibration  data  recorded  by  the  MSS, 
but  since  the  water  attenuation  algorithm  does  not  require  absolute 
radiances  this  calibration  was  not  attempted. 


TABLE  1 

SPECTRAL  BANDS  RECORDED  BY  M-8  SCANNER  SYSTEM 


ignat ion 

Wavelength  (urn) 

C9 

.43-. 52 

C8 

. 50-.54 

C7 

.52-. 57 

C6 

.55-. 60 

C5 

k-n 

00 

1 

O' 

4> 

C4 

.62-. 70 

C2 

.67-. 94 

A1 

thermal  IR 

The  first  step  in  the  analysis  procedure  is  to  determine  the  deep¬ 
water  signals.  This  is  accomplished  by  averaging  the  signals  over  a 
portion  of  the  scene  within  which  the  contribution  due  to  bottom 
reflectance  is  negligible.  For  the  scene  under  consideration  such  an 
area  was  easily  located  off  the  edge  of  the  Bank  where  the  water  depths 
are  on  the  order  of  several  hundred  meters. 

Having  determined  the  deep-water  signals  the  X^  variables  were 
computed  and  scatter  plots  similar  to  those  shown  in  Figure  1  were 
generated  for  each  band  pair.  Scatter  plots  were  generated  for  several 
areas  in  an  attempt  to  locate  an  area  of  homogeneous  bottom,  as 
indicated  by  the  linear  distribution  shown  in  Figure  lb,  which  was 
obtained  from  an  area  just  west  of  North  Cat  Cay.  Ratios  of  the  water 
attenuation  coefficients  were  then  calculated  as  outlined  in  the 
preceding  section.  The  results  are  shown  in  Table  2.  Unfortunately, 
no  in  s_itu  water  optical  measurements  were  made  at  the  time  of  the 
overflight,  but  the  ratios  shown  in  Table  2  are  within  the  range 
expected  for  oceanic  water. 


TABLE  2 

RATIOS  OF  WATER  ATTENUATION  COEFFICIENTS 
DETERMINED  FROM  AIRCRAFT  MSS  DATA 


Band  Pair 
C9/C4 
C8/C4 
C7/C4 
C6/C4 
C5/C4 


Ratio  of  At t enuation  Coefficients 
0.28 
0.29 
0.33 
0.37 
0.64 


Analysis  of  tile  Landsat  data  set  for  frame  2993-14335  proceeded 
in  a  similar  manner  to  the  aircraft  data  set.  Deep  water  signals  were 
calculated  for  an  area  off  the  edge  of  the  Bank  west  of  North  Cat  Cay, 
and  scatter  plots  of  (11SS  4)  versus  (MSS  5)  were  generated  for 
several  areas.  Due  to  the  coarser  spatial  resolution,  areas  of 
uniform  bottom  type  are  somewhat  more  difficult  to  find,  and  the  set 
of  useable  data  points  is  smaller.  Figure  3  shows  a  scatter  plot  for 
an  area  south  of  North  Cat  Cay  which  appears  to  be  relatively  homo¬ 
geneous.  The  ratio  of  attenuation  coefficients  (K4/K5)  obtained  from 
this  set  of  points  is  0.24. 

The  results  of  the  water  optical  parameter  determination  may  be 
evaluated  by  a  comparison  with  field  measurements  of  the  irradiance 
attenuation  coefficient.  Unfortunately,  no  exactly  simultaneous  field 
data  are  available.  However,  as  mentioned  earlier,  two  measurements 
were  made  near  the  test  site  within  a  week  of  the  Landsat  overpass, 
which  indicate  at  least  the  range  of  attenuation  coefficients  occurring 
there.  The  first  measurement  was  made  on  October  4,  1977  just  south  of 
North  Cat  Cay  (Station  C-5 ,  25°  32.2'N,  70°  16.9'W)  while  the  second 
measurement  was  made  on  the  following  day  west  of  Bimini  (Station  D-7, 
25°  43.9'N,  79°  13.2'W).  The  irradiance  measurement  in  the  two 
wavelength  bands  are  shown  in  Figure  4,  and  the  attenuation  coef¬ 
ficients  derived  from  the  measurements  are  shown  in  Table  3.  The 
higher  attenuation  at  0-7  is  thought  to  be  due  to  resuspension  of 
sediments  during  high  wind  conditions  occurring  on  the  4th  and  5th  of 
October.  Winds  were  calm  prior  to  the  Landsat  overpass  and  the  atten¬ 
uation  coefficients  probably  returned  to  the  values  measured  at  C-5. 

The  ratio  of  attenuation  coefficients  extracted  from  the  Landsat  data 
was  0.24,  which  is  within  ten  percent  of  the  measured  value  at  C-5. 

The  ratio  of  attenuation  coefficients  obtained  from  the  aircraft  data 
for  a  similar  pair  of  wavelength  bands  (C7  and  C4)  was  0.33,  which  is 
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F igure  3 


X2(MSS  5) 


Scatter  plot  of  X.  values  for  Landsat  bands  MSS4  (.5-. 6  Urn) 

l 

and  MSS5  (.6-. 7  Mm)  over  sand  bottom  south  of  North  Cat  Cay 
is  log  of  bottom-reflected  signal. 


Sea  Celi/Deck  Cell 


Depth  (m) 


4.  Underwater  irradiances  measured  with  filtered  photometer 
at  two  locations  on  the  Great  Bahama  Bank.  Green  filter 
corresponds  to  Landsat  band  MSS4  and  red  filter  corresponds 
to  MSS5.  Underwater  irradiances  (measured  by  "sea  cell") 
were  normalized  by  above-surface  irradiances  ("deck  cell")  in 
order  to  compensate  for  illumination  changes. 


closer  to  the  value  measured  at  Station  D-7.  Since  the  aircraft  data 
was  collected  some  10  months  later,  not  too  much  significance  can  be 
placed  on  this  comparison,  but  the  attenuation  parameters  obtained 
from  the  aircraft  data  analysis  an*  apparently  within  reasonable 
1  i  m  i  t  s . 

TABLE  3 


I RRADIANCE  ATTENUATION  COEFFICIENTS 
MEASURED  BY  PHOTOMETER  NEAR  TEST  SITE 

Creen  Filter  Red  Filter 


St  at  ion 

Date 

Locat  Lon 

(MSS4) 

(MSS5) 

Rat  lo 

C-5 

10/4/77 

25° 

32. 3' N,  79° 

16 . 9  ’  W 

.080  rtf1 

.352  m"1 

.  227 

D-7 

10/5/77 

25° 

43.9'N,  79° 

18.2’W 

.134  m"1 

.403  m  _1 

.  328 

2.2  ACTIVE  METHODS 

In  addition  to  the  collection  of  passive  data  (reflected  sun¬ 
light),  the  ERIM  M-8  scanner  is  capable  of  active  operation  using  a 
pulsed-laser  source.  The  characteristics  of  this  device  have  been 
described  by  Hasell  £t  al  [3].  Basically,  a  narrow  pulse  of  light  is 
transmitted  downward  and  the  time  history  of  the  reflected  pulse  is 
recorded  using  a  transient  digitizer.  The  transmitted  pulse  is  timed 
and  directed  to  correspond  to  a  known  pixel  location  in  the  passive 
data.  Over  shallow  water,  the  returned  pulse  has  recognizable  compo¬ 
nents  due  to  reflection  from  the  water  surface  and  from  the  ocean  bot¬ 
tom  (c.f.  Figure  5).  Thus,  the  difference  in  the  time  of  arrival  of 
these  pulses  can  be  measured  and  used  to  accurately  calculate  the 
water  depth. 

Having  a  set  of  depth  measurements  at  known  locations  in  the  pas¬ 
sive  data,  a  regression  analysis  can  be  performed  between  the  depths 
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Figure  5.  Plot  of  amplitude  of  reflected  laser  pulse  versus 

time  (horizontal  axis)  over  shallow  water  near  North 
Cat  Cay.  Incident  pulse  was  emitted  vertically  downward 
from  aircraft.  Surface-reflected  pulse  is  on  the  left, 
bottom-reflected  pulse  is  to  the  right.  Time  scale  is 
20  ns/division. 


ami  tin.'  passive  signals  in  order  to  extract  absolute  water  attenuation 
coo  It icients  tor  the  passive  bands.  As  in  the  case  of  the  passive 
method  described  in  section  2.1,  this  procedure  required  a  set  of 
points  over  a  fairly  uniform  bottom,  which  can  be  identified  by  an 
examination  ot  the  internal  consistency  of  the  data  aided  by  photo 
interpretation  or  field  observations  of  the  scene. 

2'his  procedure  was  tested  using  data  from  a  second  aircraft  run 
over  the  North  Cat  Cay  test  site  at  an  altitude  of  300  meters  with  the 
pulsed  laser  in  operation.  Water  depths  were  obtained  from  the  pulsed 
laser  data  over  a  set  of  30  points  identified  as  sand  south  of  North 
Cat  Cay.  I'lie  passive  data  values  at  the  corresponding  locations  were 
then  extracted  and  the  X^  variables  (c.f.  equation  1)  calculated  for 
each  hand.  Plots  of  these  variables  versus  depth  for  three  of  the 
passive  bands  (C3,  C6 ,  and  C4)  are  shown  in  Figures  b  (a)  -  (c) . 
Regression  analyses  were  carried  out  for  each  band  to  vield  a  set  of 
equat ions 

X.  =  a.-b.Z  (5) 

i  it 

where  i  is  the  band  number.  On  the  basis  of  the  reflectance  models 
discussed  in  earlier  reports  [1,21,  the  slope  of  the  regression  line 
(b.)  may  be  taken  as  approximately  twice  the  irradiance  attenuation 
coefficient  (K.).  The  values  of  K.  obtained  from  this  analysis  are 

it 

listed  in  Table  4  and  are  plotted  versus  wavelength  in  Figure  7.  Also 
included  on  this  plot  are  the  attenuation  coefficients  tabulated  by 
Jerlov  [4]  for  water  types  II  and  III,  and  the  measured  values  discus¬ 
sed  in  section  2.1.  Agreement  among  these  sets  of  values  is  quite 
good,  although  the  increase  in  K  with  wavelength  is  somewhat  smaller 
for  the  remote  sensing  values  than  for  the  in  s_it_u  measurements  or  the 
Jerlov  values.  The  ratios  of  the  K  values  shown  in  Table  4  are  also 
quite  consistent  with  those  obtained  passively  (c.f.  Table  2) 
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Figure  6a.  Passive  signals  in  band  C8  versus  depth.  Vertical  axis 
is  log  of  bottom-reflected  signal.  Line  indicates 
least-squares  fit.  Scatter  of  measured  values  is  due  to 
differences  in  bottom  reflectance. 


DEPTH'  M  ' 


Passive  signals  in  band  C6  versus  depth.  Vertical  axis  is 
log  of  bottom-reflected  signal.  Line  indicates  least-squares 

fit.  Scatter  of  measured  values  is  due  to  differences  in 
bottom  reflectances. 
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Figure  6c.  Passive  signals  in  band  C4  versus  depth.  Vertical  axis  is 
log  of  bottom-reflected  signal.  Line  indicates  least- 
squares  fit.  Scatter  of  measured  values  is  due  to 
differences  in  bottom  reflectance. 
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TABLE  4.  IRRADIANCE  ATTENUATION  COEFFICIENTS  FOR  WATER  OBTAINED 
FROM  M-8  SCANNER  DATA  COLLECTED  NEAR  NORTH  CAT  CAY  ON  8/12/78 


Band 


Wavelength  (nm) 


K.  (m  1) 
1 


Regression 
Coefficient  (r) 


C9 

.48-. 52 

.09  3 

0.77 

C8 

.50-. 54 

.  108 

0.81 

C7 

.52-. 57 

.114 

0.82 

C6 

.55-. 60 

.  146 

0.90 

C5 

.58-. 64 

.203 

0.92 

C4 

.62-. 70 

.  335 

0.95 

19 


_  0.4 


a 


o 

c 

o 


0.02 


0.01 


«  Low  aircraft  data  (8/12/78).  K  values  obtained 
by  regression  of  passive  signals  vs.  depths 
obtained  from  pulsed  laser. 

a  In  situ  photometer  meastjrement  (10/4/77) 
near  N.  Cat  Cay 

o  In  situ  photometer  measurement  (10/5/77) 
near  Bimini 


0.5 


0.6 

Wavelength  (gm) 


0.7 


Figure  7.  Comparison  of  irradlance  attenuation  coefficients 

obtained  from  M-3  scanner  data  with  values  measured 

[4  1 

in  situ  and  coefficients  tabulated  by  Jerlov 
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51 


The  second  type  of  information  which  can  be  obtained  from  the 
pulsed  laser  data  is  the  effective  attenuation  coefficient  for  the 
laser  beam.  This  coefficient  may  be  obtained  by  a  regression  analysis 
of  the  amplitude  of  the  bottom-ref lected  pulse  against  the  depth.  A 
plot  of  the  natural  logarithm  of  this  amplitude  versus  depth  is  shown 
in  Figure  8  for  a  set  of  sand  points  south  of  North  Cat  Cay.  The  slope 
of  the  regression  line  for  this  data  is  -0.618,  with  an  r-value  of  0.94. 
Thus,  the  effective  attenuation  coefficient  is  0.309  m  ^  at  the  laser 
wavelength  of  0.52  pm.  No  i_n  s_itjp  measurements  of  the  beam  attenuation 
coefficients  were  made,  but  it  is  interesting  to  note  that  the  empir¬ 
ical  relationship  reported  by  Shannon  [5]  between  the  beam  attenuation 
coefficient  (a)  and  the  irradiance  attenuation  coefficient  (K)  at 
0.535  pm, 

K  a  0. 2a+  0.04  m-1  '  (6) 

holds  quite  accurately  for  the  laser  attenuation  coefficient  and  the 
irradiance  attenuation  coefficient  obtained  from  the  passive  band 
CG  (0.05-0.54  pm).  Substituting  the  value  0.309  m  ^  for  a  in  equa¬ 
tion  (6)  yields  the  value  0.102  m  ^  for  1C,  which  compares  quite  closely 
to  the  value  of  0.108  m  ^  in  table  4  for  band  C8. 
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Figure  8.  Amplitude  of  bottom-reflected  laser  pulse  versus 
depth.  Line  indicates  least-squares  fit  to  data. 
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BOTTOM  FEATURES 

3.1  DATA  ANALYSIS 

The  reflectance  of  shallow  water  areas  to  solar  illumination  is  a 
function  of  the  water  depth,  the  water  optical  properties,  and  the 
bottom  reflectance.  Assuming  the  water  optical  properties  to  be 
uniform  over  a  given  scene  area,  the  signals  recorded  by  a  multispec- 
tral  scanner  system  may  be  combined  to  obtain  information  on  the  bottom 
reflectance  without  knowledge  of  the  water  depth.  The  method  of 
extracting  this  information  has  been  described  in  previous  reports 
[1,2]  and  has  been  evaluated  for  a  test  site  in  relatively  turbid 
water  [6|  using  low  altitude  aircraft  MSS  data.  An  evaluation  of  the 
same  algorithm  in  clear  oceanic  water  is  described  in  this  report  using 
higher  altitude  aircraft  and  Landsat  data.  The  analysis  of  the  MSS 
data  is  described  in  this  section,  and  an  evaluation  of  the  results  by 
comparison  with  field  observations  is  presented  in  the  following 
sect  ion. 


The  basis  for  the  bottom  reflectance  algorithm  is  the  behavior  of 
the  variables  discussed  in  Section  2.1  and  illustrated  in  Figure  1: 
namely,  that  for  a  given  bottom  reflectance  the  data  points  fall  along 
a  straight  line  whose  slope  is  K^/K.,  while  changes  in  bottom  cause  a 
displacement  perpendicular  to  this  line.  We  may  therefore  define  the 
variable 


K.  ln(L.  -  L  .)  -  K.  ln(L  -  L  ) 

Y  .1 _ *  81  1  ]  8]  (7) 

as  a  measure  of  this  displacement.  Note  that  equation  (7)  can  be 
re-written  so  that  it  involves  only  the  ratio  of  attenuation  coef¬ 
ficients  Ki/Kj,  which  can  be  obtained  by  the  procedure  described  in 
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section  2.1.  Under  the  assumption  that  the  bottom-reflected  radiance 

(L.-L  )  is  proportional  to  the  bottom  reflectance  and  exponentially 

1  s  1 

dependent  on  the  water  depth,  the  variable  is  independent  of  the 
water  depth  and  is  related  to  the  bottom  reflectance  as  follows: 


K  In  r  -  K  In  r 
j  t  i  .1 


(8) 


whe  re 
f  i  xe  d 


r.  is  the  bottom  reflectance  in  band  i  and  Y,  is  a  constant 
i  10 

illumination  and  atmospheric  conditions. 


for 


This  procedure  can  be  implemented  as  an  MSS  data  processing 
algorithm  by  calculating  the  variable  Y^  at  each  point  in  the  scene 
and  using  this  variable  as  a  depth-invariant  index  of  the  bottom  type. 

If  only  two  wavelength  bands  are  involved,  this  variable  would  repre¬ 
sent  the  end-point  of  the  processing  and  could  be  used,  for  example, 
to  generate  a  film  image  on  which  the  density  is  proportional  to  the 
Y.  value.  This  image  could  then  be  interpreted  as  a  map  of  the  bottom 
reflectance,  in  which  the  effects  of  water  depth  variations  have  been 
removed.  Alternatively,  this  index  could  be  calculated  for  a  number  of 
band  pairs  and  the  resulting  set  of  values  could  be  used  as  inputs 
to  a  multispectral  pattern  recognition  or  classification  routine  in 
order  to  categorize  the  bottom  into  a  set  of  discrete  classes.  In 
that  case,  the  algorithm  could  be  viewed  as  a  preprocessing  step  to 
remove  the  effects  of  water  depth  variations,  analogous  to  existing 
routines  for  removing  atmospheric  effects  or  scan  angle  variations. 

The  aircraft  data  set  set  described  in  Section  2.1  was  processed 
using  this  algorithm  in  order  to  evaluate  the  performance  of  the 
algorithm  under  clear-water  conditions.  After  preparing  the  data 
set  as  described  in  section  2.1  and  determining  the  ratios  of  attenuation 
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coefficients  for  each  band  pair,  the  next  step  is  to  choose  a  set  of 
wavelength  bands  for  processing. 

The  choice  of  wavelength  bands  to  be  used  in  the  bottom  reflect¬ 
ance  processing  is  determined  by  a  trade-off  between  depth  of  penetra¬ 
tion  and  sensitivity  to  reflectance  changes.  For  a  small  (Ar)  in  the 
bottom  reflectance  (r),  assuming  the  fractional  change  to  be  the  same 
in  both  bands,  the  change  in  the  bottom  index  is 


(9) 


Thus,  the  sensitivity  to  changes  in  bottom  reflectance  is  largest  for 
bands  having  the  greatest  difference  in  attenuation  coefficients.  On 
the  other  hand,  this  difference  cannot  be  indefinitely  increased 
without  simultaneously  decreasing  the  depth  of  penetration  (approxi¬ 
mately  1./K.,  if  K.  is  the  larger  of  the  attenuation  coefficients). 
Obviously  one  band  should  be  chosen  to  minimize  K. ,  but  the  other  band 
must  be  chosen  to  achieve  a  compromise  between  sensitivity  and  pene¬ 
tration  depth.  In  the  data  set  under  consideration,  data  quality  is 
also  a  consideration,  since  the  gain  settings  of  some  bands  were  either 
higher  or  lower  than  optimal.  For  this  reason,  band  C7  was  chosen 
instead  of  C8  or  C9  on  the  basis  of  data  quality,  although  these  bands 
had  a  slightly  smaller  attenuation,  and  band  C5  was  chosen  to  balance 
sensitivity  against  penetration  depth. 


Having  selected  the  wavelength  bands  and  determined  the  ratio  of 
water  attenuation  coefficients  for  these  bands  (c.f.  Table  2),  the 
bottom-type  index  was  calculated  point-by-point  for  the  scene,  and  an 
output  image  was  generated  from  this  index.  This  image  is  shown  in 
Figure  9,  along  with  the  single-band  images  for  bands  C7  and  C5. 
Exposed  land  areas  were  edited  out  using  the  thermal  band,  and  appear 


25 


collected  on  12  August  1978  at  3.3  km  altitude. 


as  white  in  these  images.  The  image  brightness  in  Figure  9c  may  be 
interpreted  as  being  approximately  proportional  to  the  bottom  reflec¬ 
tance  (neglecting  spectral  variations  in  the  reflectance),  whereas 
Figures  9a  and  9b  are  influenced  by  water  depth  as  well  as  bottom 
reflectance.  Thus,  the  darker  areas  on  the  left  side  of  the  single¬ 
band  images  appear  as  a  brighter  tone  on  the  bottom  image,  since  these 
areas  are  characterized  by  a  large  water  depth  but  a  fairly  highly 
reflecting  bottom.  The  dark  areas  near  the  islands,  on  the  other  hand, 
appear  dark  on  all  the  images  since  they  correspond  to  vegetated  areas 
in  shallow  water.  A  further  evaluation  of  these  results  is  presented 
in  section  3.2  of  this  report. 

Next,  the  Landsat  data  set  (frame  2293-14385)  was  processed  in 
order  to  evaluate  the  utility  of  the  bottom  reflectance  algorithm  as 
applied  to  data  collected  from  satellite  altitudes,  at  a  coarser  spatial 
resolution.  Using  the  deep-water  signals  and  the  ratio  of  attenuation 
coefficients  determined  from  the  analysis  described  in  section  2.1  the 
bottom-type  index  (Y  )  was  generated  point-by-point  for  the  portion  of 
the  frame  covering  the  North  Cat  Cay  test  site.  The  image  generated 
from  this  index  is  shown  in  Figure  10,  along  with  the  raw  data  images 
for  MSS  4  and  MSS  5  .  Land  was  edited  out  from  these  images  using  the 
signal  in  MSS  7  (0.8-1. 1  pm)  to  discriminate  between  land  and  water. 

A  small  geometric  correction  was  also  applied  to  this  data  to  correct 
for  image  skew  due  to  earth  rotation  effects. 

Making  allowances  for  the  difference  in  spatial  resolution,  the 
Landsat  bottom  map  (Figure  10c)  is  similar  to  the  aircraft  bottom  map 
(Figure  9c)  in  the  shallow-water  areas  near  the  islands.  In  deeper 
water,  on  the  left  side  of  the  scene,  the  Landsat  image  is  generally 
noisier  and  has  a  higher  porportion  of  dark  bottom  classifications. 

This  is  probably  due  to  the  relatively  higher  water  attenuation  in 
MSS5  than  in  band  C5  of  the  aircraft  scanner,  and  the  lower  signal-to- 
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noise  ratio  of  the  Landsat  scanner.  There  is  also  a  dark  spot  in  the 
Landsat  image  south  of  North  Cat  Cay,  which  is  caused  by  the  shadow  of 
a  small  cloud  appearing  in  the  lower  right  corner  of  the  image. 

Further  evaluation  is  presented  in  the  following  section. 


3.2  EVALUATION  OF  RESULTS 

The  results  of  the  bottom  reflectance  mapping  algorithm  were 
evaluated  by  direct  observations  of  the  bottom  conditions  at  nine 
stations  within  the  test  site.  The  locations  of  the  stations  are 
shown  in  Figure  2,  and  a  summary  of  the  observations  is  contained  in 
Table  5.  The  bottom  reflectances  indicated  in  Table  5  were  obtained 
by  analyzing  underwater  photographs  of  the  bottom,  examples  of  which 
are  shown  in  Figure  11.  Each  photograph  included  a  calibrated  gray 
scale  having  reflectances  of  1.8  percent,  8.6  percent,  and  28.7  percent 
for  each  panel.  The  photographs  were  scanned  with  a  microdensitometer 
to  determine  the  film  density  for  each  gray  panel  and  for  the  bottom. 
The  gray  panel  readings  were  used  to  obtain  a  relationship  between 
fil  rensity  and  panel  reflectance,  which  was  then  used  to  infer  the 
bottom  reflectance.  Adequate  sampling  of  the  bottom  reflectance  is  a 
problem,  since  the  field  of  view  of  the  underwater  camera  is  much 
smaller  than  the  resolution  of  the  remote  sensing  data.  In  order  to 
alleviate  this  problem,  attempts  were  made  to  locate  the  stations  in 
relatively  homogeneous  areas,  with  the  exception  of  stations  A-l  and 
B-3  which  were  located  at  the  boundaries  of  vegetated  areas. 

After  processing  the  aircraft  and  Landsat  data,  the  coordinates 
of  each  station  were  determined  by  overlaying  a  map  (Figure  2)  on  the 
processed  image.  The  bottom-type  indices  (the  Y-values)  were  extracted 
from  the  digital  data  at  each  station,  and  are  shown  in  Table  6  and 
Figure  12.  For  the  aircraft  data,  the  Y-value  predicts  the  bottom 
reflectance  with  a  standard  error  of  1.8  reflectance  units.  For  the 
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SUMMARY  OF 

TABLE  5 

BOTTOM  OBSERVATIONS 

AT  TEST  LOCATIONS 

at  ion 

Depth  (m) 

Reflectance  (%) 

Description 

A- 1 

3 

16.0 

boundary  of  vegetated  area 

A- 2 

3 

5.0 

Thalassia  bed 

A- 3 

3 

32.0 

white  carbonate  sand 

B-l 

5 

16.5 

hard,  non-vegetated  bottom 

B-2 

4 

30.0 

white  carbonate  sand 

B-3 

3 

12.5 

boundary  of  vegetated  area 

B-4 

4 

5.0 

Thalassia  bed 

C-l 

3 

35.0 

white  carbonate  sand 

C-2 

4 

3.0 

Thalassia  bed 

30 


TABLE  6 


BOTTOM-TYPE  INDICES  CALCULATED  FROM 


A I RCRAFT 

AND 

LANDSAT  DATA 

AT  TEST 

LOCATIONS 

Aire ra ft 

Data 

i 

Landsat  Data 

Stat ion 

C7 

C5 

Y1 

MSS  4 

MSS  5 

h 

A- 1 

L6L 

162 

2.  1 1 

96 

35 

3.26 

A- 2 

77 

77 

1.53 

59 

23 

2.48 

A- 3 

248 

237 

2.40 

109 

27 

3.60 

B-l 

120 

89 

2.  14 

87 

22 

3.41 

B-2 

L92 

150 

2.37 
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2 . 22 
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Values  tabulated  under  C7,  C5,  MSS4,  and  MSS5  are  signal  values 
(in  digital  counts).  is  bottom-type  index  calculated  from 

equation  (7). 
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Figure  12.  Measured  Bottom  Reflectance  Versus  Value; 
from  Aircraft-  and  Landsat  Data. 
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Landsat  data  the  standard  error  is  3.6  reflectance  units  for  the  nine 
measurements.  It  may  he  noted  that  over  the  range  of  depths  for  which 
bottom  reflectance  measurements  were  made,  the  single-band  radiances 
are  fairly  good  predictors  of  the  bottom  reflectance.  This  is  not 
expected  to  be  true  over  a  wider  range  of  depths,  however,  since  water 
attenuation  effects  would  cause  a  larger  variation  in  the  observed 
radiance  than  bottom  reflectance  effects.  To  test  the  expectation,  an 
area  further  offshore  was  selected  (see  box  in  Figure  2)  which  has  a 
depth  of  approximately  15  meters  and  a  bottom  composed  of  white  coral 
and  sand,  according  to  DMA  chart  No.  26324.  Color  aerial  photography 
seems  to  indicate  a  bottom  similar  to  that  at  station  B-l ,  although 
the  difference  in  depth  makes  the  bottom  type  somewhat  difficult  to 
judge.  Single-band  data  values  and  Y-values  were  extracted  from  both 
the  aircraft  and  Landsat  data  sets  for  this  area,  and  are  shown  in 
Table  6.  Using  a  regression  fit  between  the  single-band  data  values 
and  the  bottom  reflectances  at  the  nine  stations,  the  bottom  reflec¬ 
tances  predicted  by  C5  and  C7  are  2.8  percent  and  6.7  percent,  respec¬ 
tively,  both  of  which  are  much  too  low.  The  bottom  reflectance 
predicted  by  the  Y-value  for  the  aircraft  data  is  21.8  percent,  which 
is  probably  quite  close  to  the  correct  value  judging  from  the  reflec¬ 
tance  at  B-l.  For  the  Landsat  data  tiie  MSS-5  radiance  is  not  signifi¬ 
cantly  above  the  deep-water  radiance,  so  the  bottom  reflectance  cannot 
be  reliably  calculated.  The  Y-value  indicated  in  Table  6  for  the 
Landsat  data  implies  a  reflectance  of  7.4  percent,  which  is  in  fact 
much  too  low. 
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4.0 

RADIATIVE  TRANSFER  MODELING 

4.1  WATER  MODELING 

A  radiative  transfer  model  consisting  of  a  modification  of  the  OSS 
approximation  [2,7]  was  used  for  the  development  and  theoretical  evalua¬ 
tion  of  the  bottom  reflectance  algorithm.  This  model  has  the  virtue  of 
computational  efficiency  and  gives  acceptably  accurate  results  for  most 
cases  of  interest  for  passive  remote  sensing,  although  it  is  known  to 
break  down  in  highly  turbid  water.  The  assumptions  of  the  model  are  that 
the  water  properties  are  vertically  homogeneous,  that  the  water  and  bot¬ 
tom  are  horizontally  uniform,  and  that  the  source  of  illumination  is 
constant  and  horizontally  uniform.  These  assumptions  are  not  overly 
restrictive  for  modeling  the  passive  case,  but  are  not  applicable  to 
the  active  (pulsed-laser)  case. 

In  anticipation  of  the  need  for  modeling  the  response  of  the  active 
scanner,  a  Monte  Carlo  simulation  model  is  being  developed.  The  advan¬ 
tage  of  this  model  is  that  it  is  "exact"  in  the  sense  that  results  can  be 
obtained  to  any  required  degree  of  accuracy  given  enough  computer  time, 
and  that  it  can  be  modified  to  simulate  complex  geometries  and  boundary 
conditions.  It  can  alsogenerate  time-dependent  solutions  and  account  for 
narrow  beam  illumination  conditions.  The  disadvantage  of  the  Monte  Carlo 
method  is  that  a  large  amount  of  computer  time  may  be  required  to  generate 
accurate  results  for  some  situations.  For  this  reason,  efforts  are  also 
being  made  to  develop  approximate  analytical  models  which  are  applicable 
to  the  active  case.  The  Monte  Carlo  program  will  be  used  to  test  the  ac¬ 
curacy  and  range  of  validity  of  these  models. 

The  method  used  in  the  Monte  Carlo  model  is  to  "follow"  a  large  num¬ 
ber  of  individual  photons  as  they  pass  through  the  medium.  The  path  of 
each  photon  is  simulated  by  choosing  the  path  lengths  and  scattering  angles 
from  the  appropriate  distribution  functions.  When  the  photon  encounters 
a  boundary,  the  appropriate  boundary  conditions  are  applied,  or  if  the 
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photon  emerges  from  tho  medium  the1  relevant  properties  (such  as  angle, 
position,  and/or  time)  of  the  emerging  photon  are  recorded.  The  number 
of  interactions  undergone  bv  the  photon  before  emerging  is  also  recorded. 
Since  the  probability  oi  scattering  (as  opposed  to  absorption)  at  each 
interaction  is  given  by  the  single-scattering  albedo,  w  ,  the  emerging 
photon  is  assigned  a  statistical  weight  (w^)'1  where  n  is  the  number  of 
interactions  undergone  before  emerging.  Thus,  the  same  set  of  statis¬ 
tics  may  be  used  to  calculate  the  reflectance  for  a  number  of  different 
values  of  m^.  Other  techniques  are  possible  to  increase  computational 
efficiency,  but  generally  at  the  cost  of  reducing  the  flexibility  of  the 
mode  1 . 

A  FORTRAN  source  listing  for  *die  basic  Monte  Carlo  model  is  given 
in  Appendix  A.  This  version  records  only  the  polar  angle  of  the  emerging 
photons  and  considers  onLy  one  boundary,  at  Z=0.  Thus,  t lie  result  of  this 
program  is  the  azimuthally  averaged  reflectance  f or  a  semi-infinite 
medium.  Modifications  for  the  pulsed  laser  operating  in  shallow  water 
would  be  to  record  also  the  position  and  time  of  the  emerging  photons 
and  to  consider  a  second  boundary  at  the  ocean  floor.  Modifications 
could  also  be  made  to  include  a  vertical  stratification  of  the  water 
p ropert ies. 

4.2  ATMOSPHERIC  MODELING 

The  primary  difference  between  the  radiative  transfer  problems  in 
water  and  in  the  atmosphere  is  the  relative  amount  of  absorption.  In 
water,  the  single-scattering  albedo  w  (defined  as  the  ratio  of 
scattering  to  total  attenuation)  is  typically  around  0.5  and  rarely 
exceeds  0.9  except  in  highly  turbid  water.  In  the  atmosphere,  on  the 
other  hand,  absorption  is  negligible  in  the  visible  region  of  the  spec¬ 
trum  and  tin.'  value  of  •.,)  ^  is  1.0  for  practical  purposes.  Because  the 
QSS  model  is  known  to  break  down  as  w  >1 ,  at  least  for  large  optical 
thickness,  a  different  model  has  been  used  in  the  past  for  atmospheric 
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i'a  l  i'ii  lat  ions .  The  double-delta  model  developed  by  Turner  [8]  is  in 
relatively  wide  use  for  this  application.  However,  a  problem  has  been 
recently  reported  19]  in  that  the  reflectance  predicted  by  this  model 
is  much  too  large  in  the  anti-solar  direction.  Although  this  problem 
was  reported  for  the  case  of  optically  thick  clouds,  it  appears  to 
exist  also  for  the  path  radiiince  in  a  relatively  thin  atmosphere.  As 
a  result,  a  re-examination  of  various  models  has  been  conducted  with 
reference  to  their  performance  in  the  conservative  scattering  use  (oi  =1). 
The  QSS  model  was  included  in  this  comparison  and  found  to  yield  results 
at  least  as  good  as  the  Turner  model  in  most  cases,  although  it  does 
break  down  for  large  optical  thicknesses.  A  new  model  was  developed 
during  the  course  of  this  study  which  is  slightly  superior  to  both  the 
QSS  and  double-delta  model  for  predicting  backseat tered  radiances, 
particularly  for  optically  thick  media. 

For  the  case  of  direct  incident  irradiance  p  E  on  a  plane-parallel 

oo 

atmosphere,  the  transfer  equation  may  be  written  [10]  as 

p(M,<t>;M;4>'  )L(T,pt  ,4>'  )dp,d(f)t 

0  -1 


to 

o 


E 
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4n 


e 


U°p(b,<l>;-1‘  A  ) 

O  O 


(10) 


where  L(i,ti,i[>)  represents  the  diffuse  radiance  which  has  arisen  from  one 
or  more  scattering  processes.  The  Neumann  solution  is  obtained  by  solving 
this  equation  without  the  integral  term  (i.e.,  setting  L=0  in  the  inte¬ 
gral)  to  obtain  a  first-order  radiance  then  substituting  the  first- 

order  radiance  in  the  integral  term  and  solving  for  the  second  order 
radiance  etc.  The  n-t*'  order  radiance  in  this  series  has  the 

physical  interpretation  of  the  radiance  arising  from  n  scatterings  of  the 
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beam.  I'll  ■  difficulty  of  this  method  is  that  analytical  solutions  are 
not  obtainable  beyond  the  first  order,  so  that  numerical  integration  is 
necessary,  and  that  the  series  converges  very  slowly  for  thick  layers  as 

i.l  -1  . 

O 

In  the  QSS  model,  tile  approximate  phase  function 


p  \4>')  = 


4  ■  1 1 F  i.S  ( p-  p  '  )  iS  (  <f> "  ) 

p(M,(j>;fir<t>") 


vni'>o 

pp'-O 


(11) 


is  substituted  for  the  actual  phase  function  p(p,<j>;|J^,40  in  equation 
(10),  resulting  in  the  set  of  equations 
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for  p  -0. 


Equation  (12)  may  be  solved  in  the  single-scattering  approximation 
(by  neglecting  the  integral  term)  to  yield 


f  -(l-o)  F)x/p  -i/p 
l.',)Cr,,„*>  -  E„[c  °  °-e  ° 


(5(p+Po)<S(i()-(J)o) 


(14) 


38 


for  ii- 0.  Substituting  this  in  equation  (13)  and  solving  for  L,  we 
obtain 


(1) 


1  . 


io'u  E  -1  '/h 
o  o  o  o 

4 ir  (b+lO  ° 


(  t  '-•[  '  ) (1/p+l/p  ) 
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05) 


for  i.i  '  0,  where  i  =  ( 1  —  . i  F)  [  and  w'  =u>  /  ( 1  —to  F).  Gordon  (71  showed  that 
—  o  o  o  o 

this  approximation  yields  a  hemispherical  reflectance  which  agrees  with 

an  exact  calculation  to  within  0.5  percent  for  w  <0.6  and  to  within 

o— 

about  12  percent  for  u>  <0.85  for  the  semi-infinite  case  (t=°°)  .  It  does 

o— 

not  predict  the  forward-scattered  radiance,  however,  and  the  errors  in 

the  back-scattered  radiance  increase  as  w  +1. 

o 

In  order  to  test  the  QSS  model  for  the  atmospheric  case  (i.e., 

=  r  <°°),  calculations  of  (0 , p , <+>)  from  equation  (15)  were  com¬ 

pared  witli  exact  results  calculated  by  Hansen  [11 ]  using  the  doubling 
method.  The  phase  function  used  for  this  comparison  is  shown  in 
Figure  13.  The  reflection  function 


ir 

R(  i  ;u,4>;u  ,<f>  )  =  ---  L (o , p , <J>)  (16) 

O  E, 

o 

is  plotted  versus  p  for  the  case  p  =1  (normal  incidence)  in  Figure  14(a), 

and  for  the  case  UQ=0.5  In  Figure  14(b).  Fairly  good  agreement  is  noted 

for  the  case  of  normal  incidence  with  T  =1,  but  large  errors  exist  for 

o 

the  other  cases.  More  accurate  results  could  probably  be  obtained  by 

doing  a  second  iteration:  i.e.,  by  substituting  equations  (14)  and  (15) 

into  the  integral  term  of  the  transfer  equation  (10)  and  solving  for 

(2 ) 

I.  (r,lJ,4>).  This  would  involve  a  numerical  integration  over  ;i  and  1  , 
however,  which  detracts  from  the  simplicity  of  the  model. 
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The  approach  taken  by  Turner  f  8]  in  the  double-delta  model  is  to 
use  the  approximate  phase  function 


p'(h,  t'lli'.'t'*)  =  4iiF  '5(p-ii')iS(i|)-'>")  +  4 11 H  iS(u+p')(S(4)-<J  '+ii) 


in  the  transfer  equation  (10)  and  to  solve  for  the  first-order  radiance 


l.(l)0,H,t>)  =  U  (i)<S(M+U,  )  +  E-(t)6(u-u  ) 6 (4>— <b  +n)  (18) 

11  I  O  O  1*0  o 

o'-  J 

Solutions  for  K  (l)  and  E+(i)  have  been  given  by  Turner  [8]  for  w  =1 
and  bv  JCaehor  et  al  [  12]  for  a)  <1 .  Substituting  equation  (18)  into  the 
integral  term  of  the  transfer  equation  (10)  yields  a  second-order  solu¬ 
tion  which  is  the  one  generally  referred  to  as  the  double-delta  model. 
Again,  further  iterations  could  be  done,  but  would  require  numerical 
integrations  and  a  more  complex  computer  code. 

A  comparison  of  the  reflection  function  for  the  double-delta  model 
with  Hansen's  exact  calculations  is  shown  in  Figures  15(a)  and  15(b). 
Figure  15(a)  shows  the  overestimate  of  the  radiance  in  the  anti-solar 
direction  (0=0°),  which  increases  with  increasing  optical  thickness. 
(This  effect  is  aiso  present  for  optical  thicknesses  less  than  one.) 
Figure  15(b)  does  not  include  the  anti-solar  direction,  but  shows  that 
the  double-delta  model  generally  underestimates  the  radiance  in  other 
d i rec  t ions. 

This  analysis  of  the  limitations  of  the  OSS  and  double-delta  model 
suggests  several  other  possible  approaches  to  obtaining  an  approximate 
solution  of  the  transfer  equation.  One  such  approach  is  to  use  a  phase 
function  of  the  form 


p'(lit,|>;ii','i>')  =  Aug  6(fi-l06(i]>-(|)')  +  1-g 


o  iU  bU 

8  -  COS 'V 


B  -  c  os'V 


(a)  eQ  =  0' 


(b)  (>0  =  60° ,  <M’0  =  0° 


Figure  15.  Comparison  of  double-delta  model  with  exact  calculations 
by  Hansen  [11].  Dashed  line  indicates  double-delta  model 
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I)  an  approximate  first- 
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in  order  to  obtain  a  first-order  solution, 
order  solution  is  obtained  in  the  form 


I.(l)(!  =  K  (  L  )<S  (n+H  )  A  <■!—=: 

-  O  i 


)  +  !  (')  +  id.,  (i) 

O  I 


(20) 


using  tile  Eddington  approximation,  this  solution  ran  then  be  substituted 
back  into  the  integral  term  of  the  transfer  equation  in  order  to  obtain 
a  second-order  solution  in  closed  form. 


Substituting  the  phase  function  (19)  into  the  transfer  equation 
one  readily  obtains 


E  (I) 


C  -(l-VC/No 
E  I  e 


(21) 


In  the  Eddington  approximation,  the  terms  1.^ ( i )  and  L^(t)  are  obtained 
by  substituting  (20)  into  the  transfer  equation  and  taking  the  first 
two  moments  of  this  equation  (i.e.,  integrating  over  fi,  then  multiplying 
by  U  and  integrating).  The  result  is  the  set  of  equations 


dl. 


=  I., 


(22) 


dbi 

dl 


.•  -i  '/p 

3  (1  -in  ')L  -  2^-  K  e 

°  °  4TT  ° 


where  !  '=  ( 1  — g)  l  and  a \'= 
o  o 

equations  has  the  solution 
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where*  C  =  - - — - iL— (26) 
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This  first-order  solution  can  then  he  substituted  back  into  the  trans¬ 
fer  equation,  and  a  second-order  solution  obtained.  The  solution  is 
somewhat  more  complicated  than  the  QSS  solution  but  is  not  untractable. 
At  1=0,  the  reflected  radiance  (p>0)  can  be  written  as 


(?) 
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A  comparison  of  this  result  with  Hansen's  exact  calculation  is  shown  in 
Figures  16(a)  and  16(b).  The  agreement  is  generally  better  than  that 
for  either  the  QSS  or  double-delta  models,  although  the  results  obviously 
are  not  exact.  A  FORTRAN  source  listing  for  the  computer  program  used 
to  calculate  these  results  is  given  in  Appendix  B. 


5.0 

CONCLUSIONS  AND  RECOMMENDATIONS 

Information  can  be  obtained  from  active  and/or  passive  mul tispectral 
scanner  data  on  both  bottom  reflectance  and  water  attenuation  parameters 
in  shallow  water,  provided  that  the  assumption  of  uniform  water  pro¬ 
perties  can  be  made.  This  information  has  potential  applications  to 
studies  of  marine  biology  and  geology,  as  well  as  to  the  remote  measure¬ 
ment  of  water  depth,  as  outlined  below.  The  water  parameters  which  can 
be  obtained  from  passive-only  data  are  the  ratios  of  the  irradiance 
attenuation  coefficients  in  the  wavelength  bands  used.  Using  active 
(lidar)  techniques  in  conjunction  with  the  passive  data,  the  absolute 
values  of  the  irradiance  attenuation  coefficients  can  be  obtained  as 
well.  These  parameters  are  determined  for  the  entire  scene,  while  the 
bottom  reflectance  is  calculated  at  each  point  in  the  scene. 

The  information  extraction  algorithms  described  in  this  report  are 
relevant  to  remote  bathymetry  studies  in  several  respects.  The  water 
parameter  information  can  be  used  to  estimate  the  maximum  depth  which  can 
be  mapped  by  remote  sensing  techniques,  and  can  also  supply  input 
parameters  needed  by  these  techniques.  The  bottom  reflectance  information 
can  likewise  be  used  to  supply  input  parameters  and  to  make  corrections 
for  bottom  color  variations  in  the  water  depth  calculations. 

It  does  not  yet  appear  to  be  possible  to  devise  a  computer  algorithm 
for  distinguishing  water  parameter  variations  from  bottom  reflectance  or 
water  depth  variations,  although  it  is  frequently  possible  for  a  trained 
photo-interpreter  to  make  this  judgement,  probably  on  the  basis  of  spatial 
or  textural  differences.  The  assumption  of  uniform  water  properties 
appears  to  be  much  more  frequently  valid  than  the  assumption  of  uniform 
bottom  reflectances.  However,  more  in  situ  water  measurements  are  needed 
to  assess  the  spatial  variability  of  the  water  attenuation  coefficient 
and  to  directly  verify  the  results  of  the  water  attenuation  extraction 
algorithms. 
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Very  littLe  direct  information  is  available  on  tlie  spatial  variation 
of  the  irradiance  attenuation  coefficient.  This  lack  of  information  is 
due  in  large  part  to  the  difficulty  of  obtaining  irradiance  attenuation 
data  in  shallow  water.  The  conventional  method  of  obtaining  the  irradiance 
attenuation  coefficient  is  to  make  a  series  of  measurements  of  the 
irradiance  at  different  depths,  using  a  submersible  photometer.  This 
method  is  difficult  to  apply  in  shallow  water  because  the  problem  of 
shadowing  the  photometer  is  more  severe  and  because  the  range  of  depths 
over  which  measurements  can  be  made  is  more  limited.  Beam  attenuation 
measurements  are  easier  to  make  in  shallow  water,  because  these  measure¬ 
ments  can  be  made  with  transmissometers  that  use  an  internal  light  Source 
and  require  a  shorter  path  length.  However,  these  measurements  are  not 
easily  related  to  the  irradiance  attenuation  coefficient,  which  is  the 
relevant  parameter  for  remote  sensing,  unless  independent  measurements  of 
the  scattering  function  are  made  simultaneously.  What  is  needed,  therefore, 
is  a  means  of  making  rapid  measurements  of  both  the  beam  attenuation 
coefficient  (■<)  and  the  volume  scattering  function  3(0).  The  irradiance 
attenuation  coefficient  can  be  obtained  from  these  measurements  by  inte¬ 
grating  the  volume  scattering  function  over  the  forward  hemisphere  and 
subtracting  this  quantity  from  the  beam  attenuation  coefficient. 

Other  areas  of  research  that  should  be  pursued  include  the  extraction 
of  water  optical  properties  using  the  volume-scattered  return  from  the 
pulsed  laser,  and  an  examination  of  surface  effects.  A  modeling  effort 
could  be  undertaken  to  determine  what  water  parameter  information  can  be 
extracted  from  the  shape  of  the  returned  laser  pulse.  The  results  could 
be  initially  tested  using  existing  pulsed  laser  data  from  the  Bahamas, 

Puerto  Rico,  or  Florida.  Definitive  testing  would  require  the  acquisition 
of  a  data  set  with  simultaneous  ini  situ  water  measurements,  however. 

Surface  reflectance  and  transmittance  effects  also  deserve  further 
examination.  These  effects  are  important  for  several  reasons.  First,  in 
the  visible  region,  they  represent  a  source  of  noise  or  interference  with  the 
subsurface-reflected  light  for  most  passive  remote  sensing  applications, 
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including  bathymetry.  A  better  understanding  of  surface  effects 
would  thus  aid  in  the  development  of  more  accurate  information 
extraction  techniques.  Secondly,  active  systems  (such  as  lidar  bathymetry 
devices)  require  that  a  recognizable  surface-reflected  signal  be 
recorded,  and  it  is  important  to  design  and  operate  these  systems  so 
that  this  condition  is  met.  Thirdly,  useful  information  on  sea  state 
and  wave  parameters  may  be  contained  in  the  surface-reflected  signals 
and  improved  techniques  are  necessary  to  extract  this  information. 


49 


pr 


1 


AITKNDIX  A:  FORTRAN  Code  for  Monte  Carlo  Mode] 


1 

2 
3 
o 

5 

6 

7 

8 
9 

10 

11 

12 

13 

10 

15 
1b 

17 

16 

19 

20 
21 
22 
23 
2a 

25 
2b 
27 

26 

29 

30 

31 

32 

33 
30 
35 
3b 

37 

38 

39 
00 
«1 
02 
03 
on 
05 
06 
07 
08 
09 

50 

51 

52 

53 
50 
55 
5b 

57 

58 

59 

60 


c  monte  carlo  radIattve  tra nsfer  program 

NAME  LIST  /InPIIT/  XO.YO,ZO,THfTA.PMI,NRMAX,NlMAX,l9,NR,P,TH,NTH, 

I  NAI.MO.NNO 

C  (Xfl,Y0,Z0)»lNlTlAL  POSITION  OF  PHOTO. 

C  (THETA.  PHI)a|NtT!*|  01  PEC  1  ION  OF  PHOTON  ( DFGREE  91 
C  NPHAxsMAXIMUM  number  OF  PHOTONS  CONSIDERED 
C  NJMAXsHaX 1  MUM  number  OF  INTERACTIONS  C0N9I DERED ( <50 ) 

C  I9*1NITJAL  SEEn  FOR  RANOOM  NUMBER  GENERATOR 
C  NRPlNOFX  OF  REFRACTION  OF  MEDIUM 
C  P(  !)*SCATTERTNG  PHASE  FUNCTION  OF  MEDIUM 
C  Th(1)*SCATTERInG  angle  (DEGREES) 

C  NTH=NUMRER  OF  SCATTERING  angle S ( < 30 ) 

C  NAlaNUMRER  OF  ANGLE  INTERVALS  FOR  COUNTING  EMERGING  PHOTONS(<50> 

C  MOdloSINGLE-ScATTFRlNG  ALBEDO 

C  NW03NUMBER  OF  v*I.UES  FOR  S  T  NGLf -SC  A  T  TER  I NG  ALBEDO(<10) 

RFAL  P(30)»TH(30)#NR 
INTEGER  BIN(50,50),NEP(30) 

RFAL  R(10),rH(1O),RF(10),THB(50).H0(10) 

DATA  XO,YO,Z0/O.»0.,0./ 

DATA  THETA, pHI/0. ,0./ 

DATA  NPMAX/1 0000/ 

OAf*  NIMAX/tO/ 

DATA  13/1/ 

DATA  NR/1.333/ 

DATA  P/<I5.77,5.39, 1 .58,  .299,.  120,  .0366.  .0156,  ,  0091  5,  .00661,  ,00601  , 
1  .00732, .00629, .01 02, .0126, .0131 » ,013#, .01 36, 13*0,/ 

DATA  TH/1.,5.,  10. .20. ,30. ,95., 60,. 75. .90,, 105, , 120 , , I  35 , , 1 50 . , 

1  165. , 170. ,|75.. 180. ,13*0,/ 

DATA  NTH/17/ 

DATA  NA I/)0/ 

DATA  M0/.1,,?,,3».R,.5,.6,,7,«8,,9,1.0/ 

DATA  NWO/9/ 

INTEGFRA2  AmS/'  •/, YES/'Y'/.Nn/'N'/ 

REALaB  OASHfS/'  - '/ 

PTsARCOSC-l.O) 

TMOP I »2 . *P 1 
DEGPAD=PI/l80, 

CALL  FCVTHBC13 
10  WRITF(6,100) 

mo  FORMAT! 'OFNTER  INPUT  PARAMETERS  (NAMELIST  FORMAT  >  » • ) 

WRITE(b, INPUT! 

RFAD(5, INPUT) 
mRITE(7, InPuT) 

C  INTEGRATE  PHASf  function 
B«PINT(P, th.nthi 
WR1TE(7,102T  fl 

102  FORMATC0BACK3CATTEBING  fraction  ■ • , F7, J) 
c  calculate  anglf  INTERVALS  and  initialize  bins  TO  ZERO 
DO  12  1*1, NaI 

IMS ( I ) *ARCOS ( SORT ( { I- I . ) /NA I ) ) /PEGRAD 
DO  12  J*l,NiMAX 
BIN(I,J)»0 
12  CONTINUE 

THfl (NA I ♦ 1 )*0 , 

NPcO 

C  NP  IS  NUMBER  OF  PHOTONS  CONSIDERED  THUS  FAR 
NLP*0 

C  NLP  IS  NUMBER  OF  PHOTONS  NOT  EMERGING  AFTER  NIMAX  INTERACTIONS 
ZLP*0. 

?0  NP*NP»l 


50 


61 

62 

6} 

64 

65 

66 

67 

68 

69 

70 

7 1 

72 
7  3 
7  4 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 
109 

no 

m 

112 

113 

114 

115 

116 

117 

118 

119 

120 


ifinp.gt.npmax)  on  to  60 
c  consider  new  photon 

N  I  30 

c  reset  initial  position  and  ancle 

XPsXO 
TPs  TO 
ZPsZO 

CTsCOS(THETa*DEGRAD) 

STsSIN(THETaoOEGRAD) 

CPsCOSIPMI  *DEGRAD ) 

SPsSINIPHI  »r>t  GRAD) 

30  IF(NI.GT.NImAx)  GO  TO  50 
N I sNf ♦ J 

c  calculate  distance  d  to  next  interaction 

RD3URANO(I$i 

Ds.ALOG(RD) 

c  calculate  position  of  photon  at  this  interaction 
x«xp*o*st»cp 
Y*TP*0*ST»SP 
Z«ZP*D»CT 

c  check  to  see  if  photon  is  still  in  medium 

IFIZ.LT.O.)  GO  TO  40 

c  consider  next  path,  reset  initial  position 
32  xp»x 

TPsT 

IP»1 

C  CAtCULATE  rotation  matrix  for  previous  direction 

A3 1 *ST*CP 
A32«ST»SP 
A33«CT 

A1  l»SORT  ( t««A3t*A3I) 

*22*A33/A l 1 
A23S.A32/A1I 
A12sa3|«A23 
a  1 3*«*31 *A2? 

c  calculate  new  photon  direction  itp,ppj  relative  to  previous  direction 

RTsURAND(IS) 

CTPsPINV(RI) 

SIP«SQRT<) ..CTP»CTP) 

RPslIRANO  ( 1$) 

PPsTMCIPIoRP 

CPP«COS(PP) 

SPP*SI N (PP ) 

C  CALCULATE  DIRECTION  IN  ORIGINAL  COORDINATE  SYSTEM 
CT*Al3*STP»CPPtA23*STP»SPPTA33«CTP 
STsSQRI (1  ,-CTaCT) 

IF (SI, £0,0.)  GO  70  30 
CPslAl I»STP.CPPtA3l*CTP)/SI 
SPs(Al2»STP»CPPAA22«STP»SPPtA32»CTP)/ST 
GO  TO  30 

c  check  to  see  if  photon  is  reflected  at  interface 

40  S7AsNR*ST 

IFISTA.LI.l,')  GO  TO  42 

CT»-CT 

|s-Z 

GO  TO  32 

C  ACCUMULATE  STATISTICS  FOR  EMERGING  PHOTONS 
42  NI sN I • 1 

ST2*STA»STA 

CT2*1.-ST2 
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r  -  - 


£ 


121 

122 

12) 

120 

125 

126 
12/ 
128 
129 
1)0 
1)1 
1)2 
D) 
1)0 
1)5 
1)6 
1)7 
1)8 
1)9 
too 
101 
102 
10) 
too 

105 

106 

107 

108 
109 
150 
1  SI 
152 
15) 
150 

155 
1 59 
157 

156 

159 

160 
161 
162 
163 
160 

165 

166 

167 

168 

169 

170 

171 

172 
17) 
170 

175 

176 

177 

178 

179 

180 


NAsNAl»CT?»l 
BlNlNA,NI)*BlN!NA,Nl  )♦! 

GO  10  20 

ACCUMUIATE  STATISTICS  FOR  “LOST*  PHOTONS 
SO  NLP*NLP»1 
ZLP*ZLPAZ 
GO  TO  20 

PRINT  OUT  STATjSTICS 
60  *RIT£(7,160) 

160  FORMATl'O  ORDER  OF', IOAi •NUMBER  OF'/'  INTERACT! On *>5X, 

I  'EMf  RGING  PHOTONS' ) 

DO  60  J*l,NtMAX 

NEP(J)sO 

UO  62  Ial.NAl 

NEP ( J )*NEP ! J)*BIN(I,J) 

62  CONTINUE 

WRITEC7.I62)  J,NEP(J) 

162  FORMAT(17,»OX,HO) 

64  CONTINUE 

zlp*zlp/nlp 

hRITE<7,160)  NIMAX,NLP,ZLP 

164  FORMAT! 'ONUmBER  OF  PHOTONS  NOT  EMERGING  AFTER' » 10,  '  INTERACTIONS  * 

i  '.no,/1  average  depth  of  these  photons  *',fio.),'  units') 

WRITE(7,166)  ( HO ( I ) , I ■ 1 , NWO j 

16b  FORMAT ('OBlDl RECTI ON At  REFLECTANCE  ( A  Z I  MUTUALLY  AVERAGED)  X  100  I* 

1  /'oangle  Interval ', iox, 1  single-scattering  albedq'/ 

2  •  (DEGREES)', 2X.10F7,)) 

HRITE<7,167)  (DASHES, 1*1, NHO) 

167  FORMAT(1X,U(*-'),IOA7) 

DO  68  1*1, Nil 

1 1 =NA I -I ♦ 1 
I2=NA I -1*2 
DO  66  K«|,NwO 
H*t, 

R(K)«0. 

DO  66  Jd.NJMAX 

W*N*WO ! K ) 

H(K)*R(K)AlOO.»W*BlN! 1 1 , J)*NAI/NPMAX/P1 
66  CONTINUE 

WRITE  (7,  168',  THB(l2),THB(tl),(R(K),K*l,NW0) 

168  F0RMAT(1X,Fb.I , •  -',FS.l,lX,10F7,)) 

68  CONTINUE 

DO  72  K*l,NwO 
W*1 . 

RH(K)*0, 

00  70  J*1,NTMAX 
w*w*mo (K ) 

RM(K)*RH(K)*100, *W*NEP(J)/NPMAX 
70  CONTINUE 
WBMAMO (K ) 

HE(K)*100.»W»NLP/NPM*X 
72  CONTINUE 

WRITE <7, 170)  (RH(K),K*1,NW0) 

170  FORMAT! 'OHEmISPHERICAL  REFLFCTANCE  X  1001'/'  ' /l Ox, 1 0F7 . )) 

WRITE (7, 172)  NIHAX, (RECK), K*l, NWO) 

172  FORMAT! 'OMAxlMUM  CONTRIBUTION  TO  HEMISPHERICAL  REFLECTANCE'/ 

|  •  FROM  PHOJONS  SCATTERED  MORE  THAN', 10,'  TIHESl'/'  ' / 1 OX , 1  OF 7 , )) 
80  WRI TE  f  6, 1  BO ) 

180  FORMAT! 'OCUnTINUE  THIS  RUN?') 

READ(5< 182)  ANS 
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ia  i 

182 

FORMAT  <A1 1 

1 82 

IF(ANS,fQ,TFS1  GO  TO  84 

181 

IF(ANS.EO.Nn)  GO  TO  10 

164 

GO  TO  80 

1 85 

84 

wRITE(6,1R4) 

186 

184 

FORMATPOENtER  NUMBER  OF  ADDITIONAL  PHUTUNS  10  Bt  CONS 1 01  RED I  ' ) 

187 

RE  AD (S, 1 86 1  NADP 

188 

186 

FORMaT(IIO) 

189 

NPMAX»NPMAX»NAOP 

190 

NPsNP-1 

191 

ZLP=ZLP*NLP 

192 

GO  TO  20 

191 

End 

1 

FUNCTION  P!nT(P.TH,NTH) 

2 

RfAL  P(I01,TH(10),MU(10) 

1 

RE  AL  *  8  A(10),B(10),S(101 

4 

DEGRaD*ARCOs(«1.0 1/180, 

5 

TWOPI *160 , *DEGR AD 

6 

DO  lo  1*1, NtH 

7 

MU( I ) *COS( Th ( I ) aDEGRAO ) 

8 

IF(P(!),GT.P(NTH)1  NF*I 

9  10  CONTINUE 

10  c  fit  forward  lobEi  p«a(Ij»(i»mu(iii»*b<i)  ,  Muti)>M»MUti»i ) 

it  nfunf-i 

12  DO  12  I»l,NFl 

11  um«ALOC(P (I  j/P(I*l))/AtOSC<l,-HU<I|l/(l,-HU< !♦»,))  ) 

10  AU)*P(I)/(|.«HU(n)**B(n 

IS  12  CONTINUE 

lo  C  CONPOTF  INTEGRAL  Sfl)  from  Mini  TO  MU»HU(I) 

17  S(l)sP(l)»(l.-HU(l))/(B(l)tl.> 

18  DO  1<I  1*2, Np 

19  S<I)s9{l-l>»<P<I>»<t,-HU(t))»P  (!•!)• 

20  14  CONTINUE 

21  c  fit  rest  of  functioni  p*a(i)fb(i)»hu  *  muci)>nu>mu(ui) 

22  NTl*NTH-l 

21  DO  lfc  I *NF , nT 1 

24  B(I)*(P(I»1)*P(I)  )/  CMU(  !♦  I  )»Mli(  I )  1 

25  ATI)=P(I 1-B(I)»MU(I1 

26  16  CONTINUE 

27  B(NTMl*8(Il 

28  A(NTH)*A(I) 

29  c  compute  sin  For  rest  of  function 

10  NF1»NF*1 

H  DO  18  I*NF 1  .  NTH 

12  AMU*(MU(I-n+MU(n)/2, 

11  DMU*MU(1«1).NU(I) 

14  3(I)*S(l«n*(An*I)tB(I«l)*AMU)*DNU 

15  18  CONTINUE 

16  ST*S<n 

IT  IF(MU(I).NE,-1.)  ST*STtCAfntB(I)»(HU(n-l,J/2,)*{NU(IjTl,) 

18  C  COMPUTE  BACKSCaTTERING  FRACTION 

19  00  20  I*}, NTH 

40  IF(HuU))  2?,  22, 20 

41  20  CONTINUE 

42  22  1*I»} 

4 1  IFU.GE.NF)  GO  10  24 

44  FS«»(I)t(A(l)>P(l))/(B(I)tl.) 

45  GO  TO  26 
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?4  FS»S(1  )♦(*(!  )*Bf  n*MOf  n/?.)*HUU) 

?6  PINT*! ,-FS/ST 
RE TURN 

C  COMPUTF  mu  CORRESPONDING  TO  INTEGRAL  OF  St 
ENTRY  PINVlM) 

3  J  =  S  t  *  3  T 
on  to  i«i,nth 
IF(Sl.LE.SCl) )  GO  TO  52 
50  CONTINUE 
GO  TO  56 

52  IFd.GT.h  GO  TO  54 
C  MU>MU(1)! 

PlNVst  ..(!  ,.Mun))*St/S(11 

RETURN 

C  MU(I-n>MU>Mll(t)l 

14  PTNVsMU(l»l)MMU(I) -MU( I»t))»(9l”8(I“l))/(S(H«9CT“l)) 
RETURN 

C  Mu<MU(NTH)| 

56  PINV8MU(n-(l.tMU(n)*(M-S(I)T/CI.-StI5) 

RETURN 
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APPENDIX  B:  FORTRAN  Code  for  Approximate  Models 


1  REAL  M,MU(l9)fMU0,R!|4),TM!!9),SP!|9) 

2  PIsARCUSt-t.O) 

3  C*»TNITIAI  I7E  Hilt  1  )  VALUES 


4 

OP  8  1*1, 19 

5 

IMII j«!T-|).S. 

6 

MU(!)*cnsitH(!)»Pi/iao.) 

7 

SP(I)*SIN(Th( T)»PI/180.)*P(MU(n,-l.,0.) 

6 

8 

continue 

9 

C»*P£AD  ASYNMFTRY  FACTllR  G 

1  0 

10 

HRITEIb, 110) 

II 

1 1  0 

FORMATT'OENTER  ASYMMETRY  F ACTOR  1 ' ) 

12 

REAOIS, | 1 2)  G 

13 

It? 

F  ORMA  T ( 3F 1 0 , 5 i 

14 

BF  =  0, 

15 

DO  12  1=2,18,2 

lb 

BF=nF*SP(I-1 )A4,»SP<I)ASP(I»| > 

17 

12 

continue 

18 

BF«BF*PI/21b. 

19 

14 

WRITE(7,ll«)  fi,RF 

?.o 

114 

format! ' 1  ASYMMETRY  FACTOR  =',F7.3,',  RACKSCATTERInG  FRACTION 

?l 

1 

F7.3,'l’) 

22 

c««read  optical  Thickness  and  angle  of  incidence 

23 

lb 

MRITE(b,llb) 

24 

lit 

FORMAT! 'OENtER  OPTICAL  THICKNESS,  ANGLE  OF  INCIDENCE,  AND’, 

25 

i 

'  A  21 MUTHf ' ) 

2b 

RFAD!5,112,FND*10)  T,TH0,PHI 

27 

MU0*CO3(TM0*PI/180.) 

28 

Kffiref7. (22i  r, two, phi, i*i,  |0) 

29 

122 

FORMAT! 'OOPTICAL  THICKNESS*' ,F7, 3, 1 ,  ANGLE  OF  INCIDENCE  *•» 

30 

1 

F 7, 3, 1  DEGpEFS,  RELATIVE  A7IMUTH 

31 

1 

1  F 7 , 3 , '  OEGpEESl ,2,0THs*,I8!F5«2,',,j,F5,2) 

32 

X0*FyP(-T/Mh«) 

33 

xr,*FxP(-T»!l  ,-Gl/MUO) 

34 

C«!l ,«Gl»T*(2,/3,-MU0)*T 1 ,-XGi 

35 

C*C/(4./3.»(l.-G)*T) 

3b 

C»»LnOP  OVER  ANGLF3 

37 

00  32  1*1,1 9 

38 

U*MII(I) 

39 

X*0. 

40 

IFIU.NE.O.)  XsE*P(-T/U) 

4  I 

Rt=Mu0*2./3.»C»(C-l,)*U 

42 

R2*CC-1,)*(| ,-Gl*T 

43 

H3*(MU0»!O»r,-MU0)FP!U,»MU0,PHI)/3,  )/(MU0l(  1  ,-G)«U) 

4  4 

H! I IsRI*! | ,.X1-R2*XIR3*(1,-X«XG) 

45 

R!I)*0,75*M||0»R!I) 

4b 

32 

CONTINUE 

47 

WRITF!7,132\  (R!H,I*I,191 

48 

132 

FORMAT! 'OITFRATFD  EDDINGTON  MOOFLI'7'0  R* ' , 1 8 (F5, 3 , ' , ' ) , F 5 , 3 ) 

49 

C * *RF peat  for  turner  MODEL 

50 

DEN=4.»(Mi|0»nF*T) 

51 

DO  34  I  *  1 , 1 9 

52 

U*MU(I) 

53 

*■1. 

54 

IFIU.NE.O.)  X*1.-EXP!-T/U) 

55 

P|*P(U,MU0, PHI-180.) 

5b 

P2*P!U,-MI|0,PH1  1 

57 

Rf I»*!BF«! T-U»X)*!P| ♦P2)tMU0«X*P2)/DEN 

58 

34 

CONTINUE 

59 

NRI TE  f 7, I  3* 1  (R(I)/I*I,|9) 

bO 

134 

FORMAT! ’OTUrNFR  MODFLl'/'O  R» • , 1 8 ! F5 . 3, • , ■ ) . F5 . 3 ) 

55 
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M  C«»RfPE*T  FOR  Q3S  HrtOFL 

62  on  ss  1*1,19 

63  UaMU(l) 

69  T  S*  T*flF 

*5  rs*i. 

b6  IFCtl.FQ.O.)  GO  TO  3b 

67  TS«l,-tXP(-TS»<l,/UT|./MUO)) 

bB  36  R(I)=P(U,..M||0,PHn*YS«M<IO/(U*MUO)/BF/0, 

69  3B  CONTINUE 

HPITf(7,138)  (R(1),T*1,|9) 

T1  138  f  l)RH»T  ( '  OQSS  MOOEtl'/'O  R*  ’  ,  I B  ( F5 . 3 ,  * ,  '  )  ,  F5 , 3) 

72  GH  TO  lb 

73  END 


1  function  p(mu,miio,phji 

2  DIMENSION  Th(2S),PT(2S) 

3  DATA  TH/O., 7. 5. IS,, 22. 5, 30., 17, 5, 95., 52. 5, 60,, 67. 3. 75., 82. 5, 90., 

«  I  97. 5. 1 05., I l ?. 5, 120., 1?7. 5, 1 3S., 192. 5, ISO., 137, 5,165., 172. 5, 180./ 

5  DATA  PT / 1 OOfl ., 20 6 ., 9 . , 3 ,, 2 ,, { ,,, 5 ,, 3 ,, 2 I ,. 07 ,. 05, 

b  1  .09,  .03,  .0<i,. os,  .07,  .|a,  ,26,  .1*.  .19,  ,18,  .30,  ,50/ 

7  REAL  MU,  MIIO  ,  MR 

8  OEgRAO«ARCOS(-1. 01/180. 

9  MR*MU«MUO»SoRT(  (1  ,-MU*MII)*(l  ,-MUO*MUO)  >  oCHS  (PHI  *Df  GRAD  ) 

10  THE  T  A«ARC08 (MR) /DEGRAD 

11  DO  10  I»2, 25 

12  IFCTHFTA.LE.TH(I))  GO  TO  12 

13  10  CONTINUE 

19  P»PT(25) 

15  t2  P«PT(I)»(PTrI-l)-PT(n)«(THETA-TH(I))/(TH(I-l)-THfn) 

16  RETURN 

17  END 
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